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ABSTRACT 

Context. Circumstellar disks are characteristic for star formation and vanish during the first few Myr of stellar evolution. During this 
time planets are believed to form in the dense midplane by growth, sedimentation and aggregation of dust. Indicators of disk evolution, 
£N| ' such as holes and gaps, can be traced in the spectral energy distribution (SED) and spatially resolved images. 

Aims. We aim to construct a self-consistent model of HH 30 by fitting all available continuum observations simultaneously. New data 
sets not available in previous studies, such as high-resolution interferometric imaging with the Plateau de Bure Interferometer (PdBI) 
(*V^ at A= 1.3mm and SED measured with IRS on the Spitzer Space Telescope in the mid-infrared, put strong constraints on predictions 

■ and are likely to provide new insights into the evolutionary state of this object. 
Methods. A parameter study based on simulated annealing was performed to find unbiased best-fit models for independent observa- 
tions made in the wavelength domain A ~ l//m . . . 4mm. The method essentially creates a Markov chain through parameter space by 
comparing predictions generated by our self-consistent continuum radiation transfer code MC3D with observations. 

Results. We present models of the edge-on circumstellar disk of HH 30 based on observations from the near-infrared to mm- 
wavelengths that suggest the presence of an inner depletion zone with ~45AU radius and a steep decline of mm opacity beyond 
>140AU. Our modeling indicates that several modes of dust evolution such as growth, settling, and radial migration are taking place 
^ ■ in this object. 

Conclusions. High-resolution observations of HH 30 at different wavelengths with next-generation observatories such as ALMA and 
JWST will enable the modeling of inhomogeneous dust properties and significantly expand our understanding of circumstellar disk 
evolution. 



Key words, circumstellar matter - dust growth - planet formation - radiative transfer - edge-on disk - stars : formation, individual : 
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1. Introduction vides constraints on dust mass, grain size distribution and total 

mass. High-resolution images in scattered light from optical to 
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(N ■ Several theories have been proposed during the last decades mid-infrared (MIR) reveal attributes of small dust grains through 

; to explain the formation of planets in protostellar disks, and wavelength-dependent opacity, and interferometric continuum 

^ ■ most hypotheses are based on growth and sedimentation of dust observations in the (sub)mm regime are used to model the spa- 

' grains, e.g. the core accretion scenario (Papaloizou & Terquem tlal dust distribution of larger grains (D'Alessio et al. 2001; 

X ! 2006; Lissauer & Stevenson 2007). Nanometer-sized dust grains Dullemond & Domimk 2004; Watson et al. 2007). Resolved im- 

H - grow through low-velocity collisions by several orders of mag- a S es at different wavelengths and good coverage of the SED are 

.9.' nitude (Beckwith et al. 2000; Dominik et al. 2007; Natta et mandatory to reduce model degeneracies (Chiang et al. 2001). 
al. 2007). During this coagulation process bigger dust grains HH 30 is a young stellar object (YSO) located in the dark 

(» 1/im) decouple from the turbulent gas flow and settle in the cloud L1551 at a distance of ~140pc in Taurus. The accret- 

mid-plane of the disk (Fromang & Papaloizou 2006). This sed- ing protostar features an edge-on indirectly illuminated flared 

imentation process is expected to form a thin embedded disk disk. In the optical and near-infrared (NIR) this disk appears as 

consisting of larger aggregates (Dubmlle et al. 1995; Carballido a nebulosity separated by an obscuring belt with wavelength- 

et al. 2006). It is still debated how nature circumvents the 'me- dependent thickness (Burrows et al. 1996, hereafter B96; Cotera 

ter size barrier' that should effectively remove boulders with di- et al. 2001, hereafter C01). Interferometric observations of 13 CO 

ameters d ~lm from this inner disk by gas drag and accretion (J = 2-1) are consistent with a gaseous disk in Keplerian ro- 

onto the protostar (Weidenschilling 1977; Brauer et al. 2008) or tation around an enclosed mass of (0.45+0.04)M G that corre- 

fragmentation due to high-speed collisions (Blum et al. 2000). sponds to a typical TTauri star with spectral class M0+1 (Pety et 

Direct observations of these boulders are impossible with cur- al. 2006, hereafter P06). Observed photometric and polarimetric 

rent instruments but indirect signs of coagulation and stratifica- variability of HH 30 shows periodicity on the scale of a few days 

tion are available. The spectral energy distribution (SED) pro- and seems to be caused by either a light-house effect, periodic 



1 



D. Madlener, S. Wolf, A. Dutrey, S. Guilloteau: The circumstellar disk of HH 30 



shadowing by infalling matter, or a close companion (Watson & 
Stapelfeldt 2007; Duran-Rojas et al. 2009). The impressive dy- 
namic bipolar jet features an undulating morphology and ballis- 
tic motion that can be explained by orbital motion or precession 
of the source due to gravitational interaction with a partner sepa- 
rated by <18AU (Anglada et al. 2007). The jet is associated with 
an entrained conical molecular outflow featuring an opening an- 
gle of ~30° (P06), whose morphology is consistent with a close 
binary (Tambovtseva et al. 2008). Interferometric imaging in the 
continuum at A = 1.3mm resolved a region of reduced bright- 
ness at the center of the system, suggesting a hole in the dust 
distribution of ~40AU radius. This inner zone may have been 
depleted of dust by tidal clearing of an unresolved binary with 
a separation of ~15AU and comparable mass, e.g. Mi = O.3M 
and M 2 = 0.15M G (Guilloteau et al. 2008, hereafter G08). In this 
case, the binary components would be cooler (r e ff~3500K) and 
younger (f~lMyr) with a spectral type of ~M3. 

The aim of our campaign is to construct a model of HH 30 
consistent with all available observations from the NIR to mm- 
wavelengths. As pointed out earlier, only simultaneous optimiza- 
tion of independent observations, i.e. the well-covered SED, spa- 
tially resolved images and interferometric visibilities, can signif- 
icantly reduce model ambiguities. An analysis of the results will 
enable us to confirm signs of disk evolution, e.g. the presence of 
an inner hole, dust growth, and settling. We will also be able to 
predict possible insights with the upcoming next generation of 
observatories. 



2. Observations 

HH 30 has been studied extensively over the last decades from 
the optical band to mm- wavelengths. Most research has been 
conducted in the optical and NIR band including photome- 
try, polarimetry, medium-resolution imaging and high-resolution 
spectroscopy with ground-based instruments (Vrba et al. 1985; 
Appenzeller et al. 2005; Anglada et al. 2007) and high-resolution 
imaging with the Hubble Space Telescope (HST) from outside 
the atmosphere (B96; Stapelfeldt et al. 1999; C01; Watson & 
Stapelfeldt 2004, hereafter W04). The SED in the MIR has 
been measured with ISOCAM/ISO (Stapelfeldt & Moneti 1999; 
Brandner et al. 2000) and refined substantially by observations 
with IRS on the Spitzer Space Telescope (SST, Furlan et al. 
2008). Mapping of continuum fluxes at (sub)mm-wavelengths, 
observations of rotational transitions of CO isotopes and HCO + , 
and high-resolution imaging with interferometers have provided 
valuable insights revealing the dynamics of the gaseous compo- 
nent, its temperature distribution and the very low spectral index 
Anm~0.4 of the dust (Reipurth et al. 1993; Padgett et al. 2001; 
Moriarty-Schieven et al. 2006; P06; G08). We will present the 
most prominent observations used for this study in the following 
paragraphs. 
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Fig. 1. SED of HH 30. The data were taken from the follow- 
ing sources: Otto bolometer on KPNO: Vrba et al. (1985), 
ISOCAM/ISO: Brandner et al. (2000), IRS/SST: Furlan et al. 
(2008), MIPSJRAC/SST: Spitzer Heritage Archive, see TableQ] 
SCUBA/JCMT: Moriarty-Schieven et al. (2006), PdBI: Pety et 
al. (2006), OVRO: Padgett et al. (2001). 



Data reduction was performed with GILDAS, a radioastronomi- 
cal software package developed at IRAM (Pety 2005). Natural 
weighting was used to produce the image presented in Fig. [2] 
The beam has a FWHM extension of 0.59"x0.32" at PA = 22°, 
yielding a linear resolution along the mid-plane of ~45AU and 
~83AU in the perpendicular direction at HH 30's distance of 
~140pc. Additional details on this observation can be found in 
G08. 



2.2. IRS/SST 



2.1. IRAM Plateau de Bure Interferometer 

High angular resolution observations of HH 30 at A = 1 .3mm 
were conducted on 2007 February 3 with the IRAM Plateau 
de Bure Interferometer (PdBI) using six antennas in the A 
configuration for a total on-source observing time of 6h using 
dual polarization and a total effective bandwidth of ~1.8GHz. 
Observations of the calibrators indicated a seeing better than 
~0.15" after calibration. The rms phase noise stayed between 
19° and 53°, the brightness rms noise was 0.185mJy/beam or 
23mK with a peak brightness of 0.43K corresponding to ~19cr. 



HH 30 was observed with the IRS spectrograph on board the 
SST on 2004 February 8 using the modules Short-Low [SL] 
from 5.2-14/im and Long-Low [LL] from 14-38/im with spectral 
resolution A/AA~9Q using the stare mode (AOR key: 3552512, 
PI: J. R. Houck). The data was reduced with SMART (Higdon 
et al. 2004), which interpolates bad pixels and subtracts the sky 
(Furlan et al. 2008). The SL2 part of the spectrum was too weak 
for direct extraction and was not published, but re-reduction of 
the raw data using an off-order sky subtraction allowed sub- 
sequent extraction and calibration of the SL2 spectrum from 
1.35pm to 8.24yum (private communication with E. Furlan). 
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Fig. 2. Interferometric image of HH 30 at A= 1.3mm published 
in G08. The center lies at a = 04 A 31 m 37\469 and 6 = 
18°12'24".22 in J2000. Every contour line indicates abrightness 
step of 2<x~0.46mK with the right-hand side maximum peaking 
at 0.43K or ~19cr, the two peaks are separated by ~100AU. 



Table 1. Photometric fluxes in the mid- and far-infrared ob- 
served with the IRAC and MIPS instruments on the SST. 



A [pan] F A [mjy] 



3.6 ±0.4 


2.86 


±0.04 


4.5 ± 0.5 


2.70 


±0.11 


5.7 ± 0.7 


2.51 


±0.06 


7.9 ± 1.5 


1.68 


±0.02 


23.7 ± 2.7 


16.9 


±0.1 


71.4 + 9.5 


300.2 ± 3.9 



2.3. IRAC/SST 

The dark cloud L1551 was imaged on 2004 October 7 with the 
IRAC instrument using the map mode in all four available chan- 
nels at 3.5yum, 4.5yum, 5.8/im and 8.0yum (AOR key: 3653120, PI: 
G. Fazio). We acquired the BCD sets from the Spitzer Heritage 
Archive (SHA) v2.0 and extracted the fluxes with MOPEX 
vl8.4.9 using the Mosaic and APEX/single-source pipeline. HH 
30 was detected in all four maps as an extended source, the fluxes 
determined by aperture photometry are listed in Table Q] 

2.4. MIPS/SST 

The surroundings of HH 30 were revisited on 2006 February 16 
using the MIPS instrument in scan mode with medium rate in 
all three available channels at 24yum, 10/im and 160jt/m (AOR 
key: 12615168, PI: G. Fazio). Congruent with the analysis of 
the IRAC data we acquired the BCD sets from the SHA and ex- 
tracted the fluxes with MOPEX vl8.4.9 using the Mosaic and 
APEX/single source pipeline. Our target appeared as a point 
source in the 24/im and 70/vm maps, unfortunately it could not 
be distinguished from the background at 160/im. The resulting 
fluxes are listed in Table Q] 
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Fig. 3. Inverse images made with NIC2 using filters F110W, 
F160W and F204M aligned with the corresponding point spread 
functions (PSF) synthesized with Tiny Tim v7 . (Krist 1995). 
The PSFs were scaled by the 5th root to enhance visual contrast, 
the images are unsealed and rotated by 32° to align the disk. 

2.5. NIC2/HST 

Our object was examined by the HST on 1997 September 29 
with the NIC2 camera of the NICMOS instrument (Thompson 
et al. 1998). The observations were performed using fixed aper- 
ture and MULTIACCUM read-out mode with the filters Fl 10W, 
F160W, F187N, F204M, and F212N (HST Proposal 7228, PI: 
E. Young). This dataset has already been employed in several 
efforts to model HH 30 (C01; Wood et al. 2002, hereafter W02; 
W04). For our study we used the three images obtained with the 
mid and wide filters Fl 10W, F160W, and F204M. The three data 
sets are offered by the Hubble Legacy Archive (HLA) in its data 
release 4 (DR4), s. Fig. [3] This release applies the latest stan- 
dard pipeline for level 2 products, consisting of calibration with 
Calnica 4.4.0 and distortion correction, cosmic ray removal and 
image combination with MultiDrizzle 3.1.0. 

3. Modeling 

This section presents all models, ingredients and tools neces- 
sary to follow our approach. At the studied scale, the dust domi- 
nates the transport of radiation and thermal structure of the disk 
(Chiang & Goldreich 1997). Using one chemical mixture and 
the simple approach of Mie scattering by spherical particles, 
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our modeling effort concentrates on the spatial arrangement and 
grain size distribution of the dust. Approximating dust grains as 
perfect spheres made from one single compound may seem too 
optimistic because grains are expected to have a fractal struc- 
ture and varying composition. Then again, as has been shown 
by Voshchinnikov (2002), size, shape, and chemical composi- 
tion cannot be separated from one another, hence we restrict our 
modeling effort to the simplest shape available and one chemical 
recipe. Another necessary simplification is to neglect asymme- 
tries and variations of light sources, because the available data 
are not tangible enough to support a sophisticated model that en- 
compasses the jet, an extended protostar with accretion shocks, 
hot and cold spots, or even multiple stars without introducing 
ambiguities in abundance. With these limitations in mind we 
proceed with the description of our model. 



3.1. The disk 



The standard model for an accreting disk is the ff-viscosity 
model introduced by Shakura & Sunyaev (1973) and we use 
in this study a flared disk based on their ideas. The following 
parametrization assumes cylindrical symmetry: 



Pdust(A, z)~R a exp 



z 

'lit 1 



(i) 



This ansatz allows flaring by varying the scaling height accord- 
ing to 



MR) = h m ( — - — 
v ; U00AU/ 



(2) 



with the scaling height /iioo at a distance of 100AU and the flar- 
ing exponent f3. This canonical model has been used succesfully 
in previous studies, e.g. the Butterfly Star (Wolf et al. 2003) and 
CB 26 (Sauter, Wolf et al. 2009), but observational constraints 
require a modification. The available mm observations of HH 30 
can be explained by a hole of ~40AU radius (G08) but the SED 
in the MIR suggests the presence of warm dust (T~200K). We 
hence modified ([]]) with the step function 



6(R) 



Rin < R < Ran 

R d u < R < Rom 
else, 



(3) 



thereby introducing an attenuated inner region spanning 
from Ri n to R aa <Rout- The density jumps at this radius by the at- 
tenuation Tj € [0, 1], effectively separating the system in an inner 
and an outer disk, see Fig. |4] The dip in mm brightness observed 
by G08 is compatible with two scenarios, either shadowing of 
the inner region by optically thick layers (e.g. Wolf et al. 2008) 
or depletion of dust by a dissipation process, e.g. photoevapora- 
tion, tidal interaction with a central binary, or planet formation. 
The low peak brightness of ~0.4K in the mm map suggests the 
latter variant to be true (G08), see Fig. |2] Nevertheless, this sim- 
ple extension of the standard model is in principle compatible 
with both settings. The surface density 



^dustCR) 



-f 



PduAR, z)dz 



(4) 



follows a power law Zdust 
physics we expect 



R p with p - a - /3. From accretion 



(5) 



-200' — 
-200 
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Fig. 4. Cut through the xz-plane of the density distribution with 
depletion zone described by Eqn. ([TJ, (|2j and (0, effectively 
showing two interlaced disks. In our study the outer disk dom- 
inates the total mass and hence mm flux, while the inner disk 
contains warm dust necessary to reproduce the flux in the MIR. 
Parameters are taken from our best-fit model D, see Table |5] 

which confines the radial and flaring exponents. However, we 
treated both parameters as independent quantities, because dust 
grains can decouple from the gas stream and are subjected to 
physical processes not included in the a-disk model. Integrating 
([TJ over the entire model space yields 



K, -R ' r 

!g[g8»/jy for „ = 1 



(6) 



for the mass ratio of inner and outer disk. The total dust mass 
is adjusted to any required value »Jd U st by multiplication of ([TJ 
with an appropriate constant. Our model hence uses seven ge- 
ometrical parameters 77, a,f3, Ri D , R att , R out and /1100 to describe 
the shape and mass distribution of the disk and one additional 
parameter md us t to fix the total dust mass. 

3.2. Heating 

Three main sources of heating are known in protostellar systems, 

1 . viscous dissipation in the disk, 

2. accretion onto the protostar, and 

3. contraction of the latter. 

We neglect heating of the disk by viscous dissipation because 
the energy production in the outer parts of a protostellar disk is 
about one order of magnitude lower than the star's contribution 
(Wolf et al. 2003). Because we are not interested in asymme- 
tries, UV and X-ray excess, or line emission, we illuminated the 
system with a single black body positioned at the origin requir- 
ing two additional parameters, the effective temperature and 
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bolometric luminosity L». Owing to obscuration by the disk, the 
illuminating source of the system cannot be observed directly, 
which leads to weak constraints on these parameters. 

3.3. Dust 

It is sufficient to describe the mean grain shape, size distribution, 
and chemical composition because only mean optical properties 
of grain ensembles are used. The following simplifications were 
employed to reduce the number of free parameters necessary to 
describe a particular dust species: 

1 . spherically shaped grains, 

2. homogeneous chemical composition, and 

3. no spatial dependence of its properties. 

Given the index of refraction, the optical attributes of a homoge- 
neous sphere can be calculated with Mie's solution to arbitrary 
precision. By averaging over a given grain size distribution and 
chemical composition, a single catalog can be deduced that sum- 
marizes all optical properties necessary to model the radiation 
transfer of a particular dust species. This numerical procedure 
partitions the interval [a mm , a max ] logarithmically and evaluates 
a weighted sum (Wolf 2003). 

Size distribution The standard ansatz for modeling a grain size 
distribution is a truncated power law 



n(a) ~ a 



-d 



(7) 



with the grain size exponent d and the limiting grain sizes a mm 
and a max (Mathis, Rumpl & Nordsieck 1977, hereafter MRN). 
From these three parameters every average property of the en- 
semble can be deduced by integrating, e.g. 



{a") = 



l-d 



a"' —a 



n+\-d 



n + d — 1 



\-d_ a \-d 

mm ^ in [a max /a min ] : n = d-\. 



(8) 



The standard grain size distribution introduced by MRN to ex- 
plain the interstellar absorption has the parameter values a m ; n = 
5nm, a max = 250nm, and d = 3.5. 

Chemical composition A homogeneous mixture of smoothed 
astrophysical silicate and graphite with a ratio of 5:3 for Si:C 
and a mean density of p gra in = 2.5gcirT 3 , as described by 
Weingartner & Draine (2001), was used to model the dust com- 
position. This composition has been employed succesfully in 
the modeling of the Butterfly Star (Wolf et al. 2003) and CB26 
(Sauter, Wolf et al. 2009). We extrapolated the complex refrac- 
tion indices between 2mm<i<4mm to cover all available obser- 
vational data. Because crystalline graphite is higly anisotropic, it 
is necessary to treat the two possible alignments independently 
using the so-called "A - I" approximation 



(9) 



This approach explains the dependency of an arbitrary optical 
property Q on the components of the dielectric tensor e\\ and e± 
parallel and perpendicular to the symmetry axis and reproduces 
the observed interstellar extinction curve (Draine & Malhotra 
1993). We follow Wolf & Voshchinnikov (2004) by calculating 
the optical properties of spherical particles with the subprogram 
miex of the MC3D package using Mie's solution and empirically 
measured indices of refraction. 



Table 2. Canonical models of previous modeling work for HH 
30 except W04. All parameters marked with the subscript / were 
fixed or taken from previous publications in the corresponding 
study. Numbers marked with > are denoted as lower limits by 
the authors. W02 split luminosity into two contributions from 
the star and viscous dissipation. The dust masses assume a mass 
ratio of 1:100 between dust and gas. 





B96 


C01 


W02 


P06 


a 


2.2 


2.367, 


2.25/ 


2.2/ 


P 


1.45/ 


1.289 


1.25 


1.25/ 


flin [AU] 


0.5/ 


0.07 


0.03 


Rout [AU] 


>250, 


200, 


200, 


145 ± 20 


fcioo [AU] 


15.5 


15/' 


17 


15.5/ 


m d u St [10" 5 M o ] 


6.0/ 


>0.67 


> 1.5 


(2.7 ± 0.4) 


L,[L ] 


1.0 


0.2... 0.9 


0.2 (+0.04) 


0.2/ 


T ef! [K] 




3000 


3500/ 


3700; (Ml) 


in 


82.5 


84/ 


84/ 


81 ±3 



4. Parameter study 

The basis of understanding a phenomenon is the ability to pre- 
dict the outcome of observations by a model. A parameter study 
consists of three subsequent steps that may be iterated until suffi- 
cient agreement between predictions and observation is reached. 
As a first step a particular model and a corresponding param- 
eter space must be specified. In the following step a parameter 
set a* is identified, which reproduces the observational data best 
with the assumed model. Finally an error bound Aa for this best- 
fit parameter set is characterized to evaluate the quality of the 
model. 

Finding the most likely parameter set is equivalent to locat- 
ing the minimum of the chi-square distribution 



z 



(Ma)-y,) 2 



(10) 



derived from N measured quantities (^,+A^,,y,+Ay,) and their 
corresponding model predictions /i,(a), e.g. the SED as pairs of 
wavelength and flux (/1,+A/l;, Fi+AFi). The error is given by 



Ayf 



dxi 



(a) ■ Ax,- 



(ID 



but the second term can be neglected if \d/Ji/dxj ■ Ax,| <k |Ay,|, 
as is true most of the time. A straightforward method to iden- 
tify the global minimum is to calculate all predictions on a fine 
grid and compare them with the observational data. This exhaus- 
tive search is not feasible in high-dimensional parameter spaces 
because it is not trivial to specify a sensible grid resolution a 
priori and only a small sample of predictions can be calculated 
on a fine grid in a practical timeframe. Therefore one can only 
reduce the likelihood to find a better parameter set in most high- 
dimensional optimizations. If discrepancies between predictions 
and observations persist, the given model should be scrutinized 
and refined. This model review often introduces additional de- 
grees of freedom and this increase in dimensionality may lead to 
complications and degeneracies during the optimization process, 
see Sect. 14.41 for details of the refinements we implemented for 
this parameter study. HH 30 has already been modeled with sim- 
ilar models and less ample datasets (i.e. B96; C01; W02; W04; 
P06), see Table|2]for canonical models of these studies. 
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4.1. Simulated annealing 

We used simulated annealing (SA) to search for the best fit a» 
(Kirkpatrick et al. 1983), a Markov Chain Monte Carlo method 
(MCMC) similar to the Metropolis-Hastings (MH) algorithm 
(Metropolis et al. 1953, Hastings 1970). This method has some 
specific advantages for high-dimensionality optimization be- 
cause no gradients need to be calculated, local minima can be 
overcome intrinsically and points near optima will be evaluated 
with a higher probability. However, no upper bound for the step 
count to reach the global optimum can be given and a reasonable 
criterion for chain abortion must be defined. 

We denote in the following the sequence of all propositions, 
i.e. rejected and accepted, by (a„) and the canonical Markov 
chain of all accepted states by (a„). The random walk through 
the parameter space P - I\ X ... X Id c R d is generated by 
starting at an arbitrary point ao e P and subsequently sampling 
steps Aa from a transition probability T. We chose a Gaussian 
distribution 



r(Aa) 



Art 



(12) 



with predefined step widths /3j for every dimension. The exact 
form is not important as long as the parameter space can be ex- 
hausted in principle and a„ + Aa e P. Optimal step sizes cannot 
be determined a priori because no information is available on 
the distribution under investigation, and they have to be adapted 
along the way. After calculating all predictions ji/,(a n+ i) at the 
proposed next location in parameter space, the acceptance prob- 
ability 



A = min < 1 , exp 



A X 2 



(13) 



'^be- 



is evaluated by calculating the difference Ax = X n +\ 
tween proposed and actual position in parameter space with eq. 
( fTOb . r is the effective temperature of the "annealing" process. A 
uniformly distributed random number u e]0, 1[ is generated and 
if m < A, the step is taken, otherwise the proposition is rejected 
and a new step Aa is generated. This decision is central to the 
MH algorithm and generates the Markov chain through phase 
space by taking every downhill step while still enabling the es- 
cape from a local valley depending on the actual temperature t. 

Most implementations of SA use a monotonic cooling sched- 
ule, but these schemes contain several parameters that must be 



Table 3. Subspace for the parameter study 



parameter mm 



*1 


0.0 


1.0 


a 


2.0 


3.0 


P 


0.5 


2.0 


Rin [AU] 


0.1 


70 


#at. [AU] 


0.1 


70 


*ou, [AU] 


100 


500 


^ioo [AU] 


5 


20 


mdmt [M Q ] 


1 • io- 6 


2 • 10~ 2 


Tat [K] 


3000 


4000 


L, [L s ] 


0.1 


1.0 


flmin l^m] 


0.001 


0.1 


Omax [M m ] 


0.125 


65536 


d 


2.5 


4.0 
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adjusted empirically by running the optimization several times. 
Because the calculation of predictions p;(a) in our study is very 
time consuming, a self -regulating algorithm is preferable. We 
therefore used a non-monotonic cooling scheme by setting the 
temperature for the n-th step to 



■xt] 



(14) 



with a fixed parameter 9, the x n of the previously accepted step, 
and the xl of the global optimum (Bohachevsky et al. 1986; 
Locatelli 2000). The latter value is typically unknown and one 
has to resort to a lower bound or sensible approximation for xl 
like the^- 2 of the current best-fit of an individual run. The param- 
eter 9 < 1 has to be fixed empirically to maintain the temperature 
in the order of typical fluctuations Ax 2 encountered during opti- 
mization. 

We defined the local acceptance ratio of the chain using the 
sequence (k„) e {0, 1} of accepted and rejected steps by 



m=il-I+ 1 



(15) 



for some constant lookback length / < n. To remain in the vicin- 
ity of the optimal acceptance ratio £~0.234 (Roberts et al. 1997), 
we implemented a control loop that adjusts the step widths fij of 
the transition density (fT2l on the fly. The algorithm analyzes the 
last I positions of the chain and calculates the moduli of rela- 
tive changes between the proposed parameter vector a„ and the 
previously accepted parameter set a„_i : 



*n,j 



(16) 



The index /' e {1, ..,D] enumerates here the parameter dimen- 
sions. By separately summing up the relative change vectors 
c„ e R D in accepted and rejected proposals two vectors g„ and 
b„ can be derived, which describes the impact of good and bad 
decisions: 



s 

,,,=,1-1+1 



(17) 



If g n > f the parameter with the smallest component in g„ is 
modified to encourage riskier proposals. In the opposite case the 
parameter with the largest component in b„ is reduced to induce 
a more conservative approach. In either case the corresponding 
stepping size [3 r is adjusted proportionally by multiplying with 
(1 - £) + g n . To avoid freezing of individual parameters during 
the run, it is recommended to define a reasonable lower bound 
for every flj, i.e. enforce /3j > 6\Ij\ with 6 <sc 1. 

4.2. Calculation and combination ofx 2 

There are three fundamentally different datasets available for our 
study: the SED, resolved images, and interferometric visibilities. 
The methods to calculate the corresponding xl differ and we will 
outline the involved procedures. 

The computation of Xsed ' s straightforward because the in- 
dependent measurements can just be entered into eq. (fTOt . 

To compare observed and simulated brightness distributions 
the synthetic images must be scaled to the pixel scale of the orig- 
inal image, convolved with the PSF of the corresponding instru- 
ment setup, and then registered with the observed image. There 
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are plenty of image registration algorithms, in our case the im- 
ages only needed translation, therefore it was sufficient to evalu- 
ate the correlation integral or match the flux maxima. 

To calculate the deviation ;t^ isib for the observed visibilities 
Vi, one has first to generate visibilities Wj from a simulated im- 
age by Fourier transformation. To achieve sufficient resolution 
in the uv plane, the image must be padded with zeroes prior to 
using the FFT as size and count of pixels of an M x W array are 
connected by 



MAxAu = NAyAv = 1 . 



(18) 



The complex visibilities W, are then interpolated from the result- 
ing array using the actual antenna tracks in the uv plane. Because 
the real and imaginary parts were measured independently with 
the PdBI, we used the expression 



2 

X visib 



z 



\%m - %(wd\ 2 + mvi) - g(w,)p 



(19) 



The center of the object does not lie exactly in the origin of the 
coordinate system during observation, equation ( fT9l must hence 
be minimized by translating the visibilities in the uv plane using 
the expression 

W' t = Wi • exp [-2ni(xui + yv,-)] (20) 

to find the best translation vector (xo,yo). 

It is often necessary to combine;^ derived from independent 
observables. One straightforward method is linear combination 
of the individual values 



2 

^total 



z 



gkX 2 k , 



(21) 



which introduces an arbitrary weighting through the coefficients 
gk- Unfortunately, the signal-to-noise ratio (S/N) or even the 
noise itself can differ strongly between observables and it is not 
clear how to choose the coefficients gk- 

Confronted with this problem we devised a staggered 
Markov chain. Instead of performing all calculations at once to 
evaluate all x\ an d directly combine them to decide to accept a 
proposed configuration or reject it, we performed an MH deci- 
sion after each individual calculation of x\- We hence needed to 
keep track of individual temperatures for every dataset that is 
updated as described in the previous paragraph. By only accept- 
ing new configurations after all individual x\ steps were com- 
pleted succesfully, the algorithm effectively optimizes all obser- 
vational data sets simultaneously without weighting them. An 
added advantage of this procedure is that there is less computing 
time wasted on uninteresting parameter sets, especially if less 
time-consuming calculations are performed first. In any case, af- 
ter performing the MCMC optimization it is necessary to choose 
the weights gk to fix the metric of^ otal by analyzing the calcu- 
lated samples. 



4.3. Confidence interval 

Localizing a best-fit parameter set a* is the first step of the fitting 
process. The second step is estimating a confidence region sur- 
rounding this optimum, which can be established by sampling 
the ^-distribution and locating the confidence boundary given 
by 



Sampling of an unknown distribution was the prime motiva- 
tion behind the MH-algorithm, therefore no major code change 
is necessary to estimate error bounds. To identify an appropriate 
A-^conf the reduced merit function 



X_ 

K 



(23) 



is often used, where K is the number of degrees of freedom 
of the model. Unfortunately, x 2 ed cannot be employed for non- 
linear models because K is unknown and can change depend- 
ing on the position in parameter space (Andrae et al. 2010). 
Additionally, the S/N can differ significantly between observ- 
ables which leaves the choice of A^ onf to the practicioner's 
discretion. After choosing a sensible Ax 2 mnV the vicinity can 
be sampled by several Markov chains starting from a, with 
a constant temperature r to probe the Boltzmann-distribution 
exp [-^ 2 /tJ. The constant temperature has to be adjusted to pre- 
vent the chains from leaving the vicinity of a, too fast, e.g. 
t < Ax\ on{ ■ To bootstrap this procedure one can inspect the find- 
ings of previous Markov chains used for SA. Otherwise, the re- 
gion has to be explored first by using an educated guess for the 
temperature. After sufficient sampling, the error margin Aa* can 
be deduced by extracting the minimal and maximal components 
of all sampled parameter vectors a,- with x] < xl on f Because 
the confidence region is normally distorted by degeneracies and 
non-linearity of the underlying model, the margin Aa* is com- 
posed of non-centered intervals containing the components of 
a*, see Table [5] The reader is kindly reminded that we cannot 
exclude with certainty that there might be a similar or better so- 
lution beyond these error intervals. This would require a global 
search of the 14-dimensional parameter space on a sufficiently 
refined grid, a task not feasible with currently available comput- 
ing resources. 



4.4. Implementation 

We started our model fitting with the high-resolution interfer- 
ometric continuum images obtained with the PdBI described in 
Sect. 12.11 The dip between the two maxima in Fig.|2]is consistent 
with an inner hole in the optically thin case (G08). Molecular 
line observations made with the PdBI show a single gaseous 
disk in Keplerian rotation (P06), excluding the possibility that 
each peak of the continuum image is actually an individual disk 
revolving around one of the components of a wide binary. In the 
optically thin case, the total dust mass m^st and flux F,\ are con- 
nected by 



'Wdust'Cl = 



F A d 2 



Ba((7 'dust)) 



(24) 



-^conf AT* A^conf ' 



(22) 



with the distance d to the object and the mean temperature (7d ust ) 
of the dust. This equation can only be solved for m^st if the dust 
opacity aq is known, hence normally only the left-hand side of 
equation (l24l can be evaluated (e.g. Reipurth et al. 1993). By 
choosing a particular dust model, the total dust mass becomes 
fixed, thereby eliminating one modeling parameter in the opti- 
cally thin case. The lack of detectable emission in the interfer- 
ometric image beyond /?>130AU is inconsistent with the exten- 
sion deduced from optical and near-infrared images, as can be 
seen by superposing observations made with the HST (G08). It 
indicates a significant change in dust emissivity or drop below 
the detection limit at this distance. 

With this information in mind, we implemented a simple 
flared disk model, using the parameters a, /?, R[ n , R out , h\oo, md us t, 
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T e g, L„ and i with the standard MRN dust composition because 
there was no need for more free parameters at this initial stage of 
the investigation. By allowing for a varying dust mass the possi- 
bility of optical depth effects was acknowledged. We conducted 
a grid search on the corresponding subset of Table |3]by divid- 
ing the nine intervals of the implemented parameters into 3 ... 5 
parts. For every point on the grid the temperature distribution and 
a raytraced image at A — 1.3mm was generated, then convolved 
using a Gaussian PSF with FWHM extension 0.59"x0.32" at 
PA = 22° as described in G08 and scaled to fit the angular res- 
olution of the deconvolved mm map. The shapes of measured 
and synthesized brightness distributions were compared by nor- 
malizing to a maximum flux of 1, registering by correlation, and 
calculating ^ hape by subtracting and squaring. This approach led 
us to the following intermediate results: 

1. The interferometric image can be explained by a flared disk 
with an inner hole of R m = (55+15)AU and an outer radius 
of flout = (130±30)AU. 

2. The vertical geometry of the system can be best reproduced 
by a thin, ring-like disk with a fiat (i.e. constant) surface den- 
sity distribution. 

As a next step, we investigated the SED of best-fit models found 
by the grid search. To exclude variability of the star and un- 
known foreground extinction from the modeling process, only 
data points with A > 5/j.m were used to calculate x\et>- Soon it 
became clear that a flared disk model with a large hole and sim- 
ple MRN dust cannot reproduce the observed fluxes and spec- 
tral indices, either in the MIR or the mm domain. To counter 
the missing flux in the MIR, we introduced a depletion zone 
parametrized by an attenuation factor and radius as described 
in Sect. 13.11 effectively adding warm dust to the inner part of 
the disk. Additionally, we allowed for a variable grain size dis- 
tribution as described in Sect. 13.31 to include larger dust grains 
in the mix to improve the spectral index and total flux, espe- 
cially at longer wavelengths. Expansion of the parameter space 
by the new parameters 77, R att , a m i„, a max , and d rendered a grid 
search ineffective and led to implementation of SA as described 
in chapter l4Tl to search for the global optimum in the parameter 
subspace defined by Table [3] 

Four optimization runs were performed with this technique 
using different subsets of the observations described in Sect. |2] 
They resulted in four best-fit parameter sets that we will refer 
to as model A, B, C, and D in the following. Only the NIR im- 
ages taken with NICMOS/HST where used in the optimization 
run for model A, see Fig. [3] Model B is based on all available 
SED observations with /i>5/vm, see Fig. Q] The same data set 
was employed to find Model C, but with the restriction 77=1 to 
implement the limiting case of a simple flared disk. Model D 
is the best fit for simultaneous optimization of SED data and 
mm image from G08, see Fig. [2] For model A, we started eight 
Markov chains to fit the image observed with NICMOS using the 
F204M filter and applying in principle the same procedure im- 
plemented for the interferometric image. We rotated the original 
image by -32°, cut out a region of 4"x4" centered on our object, 
re-sized to 400 x 400 pixels and smoothed with bilinear interpo- 
lation to reduce artifacts from the rotation, see Fig. [3] Synthetic 
scattered light images were generated, then convolved with the 
corresponding Tiny Tim PSF, and finally registered and com- 
pared to the observed image. To eliminate the unknown fore- 
ground extinction and variation in luminosity, we set the error 
and peak of simulated and reference image to unity for calculat- 
ing ^p 2 04M' hence optimizing only the shape. Some fair matches 



for the F204M filter observation were found after a few hun- 
dred steps. We then calculated the corresponding images at the 
F110W and F160W filter wavelengths for the best 50 images 
found and searched for quantitative matches in all three filters 
by dereddening the simulated fluxes to match the observation in 
addition to the already fitted shape. The parameter set and^p 204M 
of the resulting model A are listed in Table [5] its optical prop- 
erties can be found in Table [4] This parameter set should not be 
understood as an overall best-fit, but just as the best model found 
by our chains. As W04 discussed extensively, the model is too 
degenerate to pinpoint a location in parameter space given scat- 
tered light images alone. Especially parameters describing the 
inner disk, i.e. R m , R an , and 77, may not be constrained at all by 
the scattered light images. 

The results of the NIR data optimization suggested that 

1 . Small dust grains with a max < 2/jm are sufficient to repro- 
duce the observed flux, optical thickness and chromaticity. 

2. Observations in scattered light are consistent with the postu- 
lated depletion zone. The specifics of the inner region do not 
change the overall appearance, as proposed by W04. 

3. Parameters influencing the vertical geometry (/zioo,/3) are 
better constrained than radial parameters (R[ n , R att , and R ont ). 

4. Contrasts in vertical cuts between the secondary maximum 
and enclosed minimum in the synthetic scatter images are 
systematically too low by a factor of ~2, see Fig. [5] 

For model B, we fitted the SED using 32 individual MCMC 
runs started from randomly selected points and propagated 
freely for ~1000 steps through the parameter subspace. Several 
Markov chains found parameter sets with low Xsed> me ^ est 
parameter set of this run is shown in Table [5] as model B. The 
low spectral index /3 mm ~0.4 indicates the presence of large dust 
grains, i.e. a max > 1mm, and a grain size exponent d < 3.5 
(Draine 2006). However, the models investigated in this param- 
eter range often lack flux and exhibit a very different spectral in- 
dex in the MIR than observed. Nevertheless, these MCMC runs 
gave us more hints: 

1. Large dust grain sizes a max > 1mm are necessary to explain 
the observed flux, especially at longer wavelengths. 

2. The flux and spectral index at mm-wavelengths can be re- 
produced with a dust exponent d~3 and a max > 8mm. 

3. No synthetic SED fits the observations at short and long 
wavelengths simultaneously. 

4. The geometry of the object is only weakly constrained by the 
SED, i.e. the model is degenerate. 

It is noteworthy that the best solution and its local confidence 
interval suggest a change in dust density by a factor of >10 at a 
distance of ~45AU, thereby dividing the system into an inner 
and outer disk, see Fig. [4] 

To ensure that this finding was not biased by introducing 
the depleted zone, we also conducted an optimization without 
this feature by fixing 77=1, which resulted in model C. We set 
up 24 Markov chains to search for solutions with this simple 
flared disk model and checked their appearance at 1.3mm. After 
~1000 steps, several Markov chains found parameter sets with 
^sed comparable to, but not better than, model B. The best pa- 
rameter set of this run is given in Table 5. The 1.3mm map gen- 
erated using these parameters is representative for the parameter 
sets for this run and is centrally peaked, showing poor agreement 
with the observed brightness profile, see Fig. [2] 

The most apparent discrepancy between models A and B is 
the different grain size distribution and thus the total dust mass. 
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Fig. 5. Contour plots and vertical cuts along the symmetry axis 
of the HST/NICMOS observations (red) and the best-fit model 
A (blue). The contours in the left column and the dotted intervals 
in the right column correspond to Icr. 

In an attempt to reconcile all approaches, we implemented the 
simultaneous optimization of SED and brightness map at A — 
1 .3mm with a staggered Markov chain as described in Sect. 14.21 
We started 24 Markov chains and let them run for -1000 steps 
to search for a unified model. The parameter set listed in Table 
[5] as model D is the best compromise. The analysis of this run 
gave us additional insights: 

1. Combination of independent datasets stall optimization, i.e. 
the Markov chains need longer to find a compromise. 

2. Staggered optimization reduces model degeneracies, e.g. the 
outer radius R out moved to values < 200AU by including the 
mm image in the optimization. 

The latter effect was expected from the results obtained by 
our initial grid search, but solutions with R out > 200AU are 
still difficult to exclude with certainty using only direct compar- 
ison of brightness maps at mm- wavelengths. To circumvent pos- 
sible artifacts introduced by the CLEAN procedure we imple- 
mented direct fitting of the visibilities in the uv plane. Analysis 
of all ~25000 synthesized mm images showed that the mea- 
sured visibilities can be best explained with R out = (135+15)AU. 
Unfortunately, these models have difficulties reproducing the 
flux in the MIR and we concluded that further optimization 
should be continued in a follow-up study when additional high- 
resolution observations in the (sub)mm domain are available that 
reveal the spatial dependence of the spectral index y6 ram . 





o 
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Fig. 6. Simulated images of model A in the NIR and the devia- 
tions from observations with NICMOS with overlayed contour 
lines spaced at Icr. As expected, the differences are lowest in 
the optimized F204M image and increase toward shorter wave- 
lengths. 

Table 4. Dust opacity k, albedo a and asymmetry factor g = 
(cos 0) for model A. Values in parentheses are from the best-fit 
model of C01. 
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1.61 


9731 (8790) 
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0.71 (0.46) 
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5. Results 

The models mentioned in the preceding section are presented in 
Table[5]for direct comparison. The parameters of previous stud- 
ies are compiled in Tableland agree well with our results. This 
is not self-evident because our parameter study was unbiased and 
employed an extended disk model with 14 parameters prone to 
degeneracies. 

The radial temperature distributions in the mid-plane of our 
models are compared in Fig. [7] with the power law deduced from 
13 CO observations by P06 



T(R) 



MlOO 



V 
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Table 5. Best-fit and deduced parameters found by optimiza- 
tion runs for the NIR images (A), SED data only (B), SED data 
only for the simple flared disk model with rj-1 (C), and mm- 
map/SED combined (D). The confidence intervals of model B 
were deduced by the procedure described in Sect. 4.3 and are 
discussed in Sect. 5.2. 
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with r 10 o = (12+1)K and q = (0.55+0.07). Vertical cuts along 
the symmetry axis of the brightness maps at A = 1.3mm are 
depicted in Fig. |8]in direct comparison with the observation by 
G08. In Fig. [9] synthetic SEDs for 5^m< A < 4mm are overlayed 
on the observed fluxes. We discuss the individual properties of 
our best-fit models in the following sections. 



5.1. Fitting of NIR observations 

Apart from small asymmetries due to uneven illumination and 
influence of the jet, model A appears to be a fair match to the 
NICMOS/HST observations, see Fig. [5] and [6] for a direct com- 
parison of simulated and measured scattered brightness distri- 
butions. However, it must be noted that we were not able to re- 
produce the contrast between the secondary maximum and the 
central minimum in the vertical cuts, as can be seen in the right 
column of Fig. [5] This phenomenon has also been reported in 
previous studies using scattered NIR images (C01 and W04), in- 
dicating a systematic deviation of our model from reality. One 
possible explanation for this discrepancy might be the cone of 
low-velocity gas driven by the central jet that has been de- 
tected in molecular lines (P06). This gas stream can raise dust 
(Tambovtseva et al. 2008) that would substantially enhance scat- 
tering of photons from the central source in this region and hence 
contrast. Again, we decided not to increase the complexity of 
the model without additional observational data supporting this 
claim because other scenarios such as disk warps, shadowing, 
and asymmetric illumination by a binary (e.g. Tambovtseva et 
al. 2006) may cause the same effect. 
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Fig. 7. Temperature in the mid-plane as predicted by our best-fit 
models and the best-fit power law of P06 deduced from 13 CO 
observations. 



The models of B96 and C01 are the counterparts of model A 
because these studies used optical and NIR observations as the 
basis for modeling. The most prominent difference is the outer 
radius R out = 455 AU, which appears to be much larger than cor- 
responding values found by our predecessors. This discrepancy 
is not as large as it may appear at first glance because the best-fit 
models of previous studies that treated R oul as a free parameter, 
i.e. B96 and P06, yielded fl out ~425AU and (420+20)AU, respec- 
tively. The flux scattered by outer parts of our models with large 
radii is much lower than lcr and is not significant compared with 
observations, especially with the relatively high surface density 
exponent p~\A of model A. Therefore, we show in Figs.|5]and 
[6]only the region with flux high enough to be meaningfully com- 
pared with observations. Our modeling results suggest that NIR 
observations alone cannot yield an upper limit for R out , therefore 
the lower values found by previous studies are lower bounds and 
our proposed value should be treated with caution. 

The ratio of dust opacities for the used filter wavelengths of 
fFiiow^Fi60w:"'F204M=l-50:l. 18:1.00 is comparable to the best- 
fit of C01 (1.37:1.29:1.00), as are most optical dust properties, 
see Table |4] 

Despite the promising appearance of model A in the NIR, its 
SED and mm map are disappointing, see Figs. [8] and [9] This 
model is missing some essential ingredient and its deviation 
from the observed temperature in the mid-plane as depicted in 
Fig.Qunderlines this fact. 



5.2. Fitting of the SED 

Model B predicts the SED quite well, see Fig.|9j and comparison 
with model A uncovers the necessary ingredient: larger grains 
with radii a > 1mm. This result is congruent with the findings 
by W02, although these authors used an exponential cutoff for 
their grain size distribution. Unfortunately, model B cannot re- 
produce the NIR observations although the overall geometry of 
models A and B are very similar. The presence of larger grains 
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in model B changes the optical properties of the dust to the 
point that this configuration is no longer a good match for the 
NIR observations. Nonetheless, these large grains are necessary 
to account for the observed fluxes and spectral indices at mm- 
wavelengths and model B is in fact a tradeoff because it is not 
able to reproduce the low spectral index /3 mm ~0.4 and the ob- 
served fluxes in the mm domain with the same confidence as in 
the MIR. Contrary to intuition, no silicate features are visible in 
the model SED although the grain ensemble contains sufficiently 
small grains. This is a geometric effect caused by shadowing of 
the inner, hot disk because tilting by a few degrees unmasks the 
hidden 10/im bright region. 

The slice through the mm map shows that the emission of 
model B at A — 1 .3mm is too extended and has insufficient peak 
brightness and contrast. Models with lower dust exponents that 
are able to match the (sub)mm flux exist, but they cannot repro- 
duce the well-sampled SED in the MIR domain. 

The radial temperature distribution in the mid-plane of 
model B follows for 7?>5AU the power law d25l > with = 19K 
and q = 0.33. The temperature varies from T~290K at the in- 
ner edge to r~10K in the outer regions with a mass average 
temperature of (7d us t) ~18K. Our model predicts slightly higher 
temperatures in the midplane than deduced from observations by 
P06, but is still a fair match (see Fig . [7J . 

Effective temperature and luminosity are consistent with 
a contracting protostar of spectral class M3+2 with R*~2R Q 
(Simon et al. 2000), but it is very difficult to constrain the prop- 
erties of the central source and a binary, as proposed by Anglada 
et al. (2007) and G08, seems a viable option. 

The surface density exponent p~1.2 is larger than expected 
for a thin accretion disk with a-viscosity, i.e. eq. (0 does not 
hold. This is no surprise because the postulated inner depletion 
zone must be generated by some additional physical process. 

The calculation of confidence intervals was performed for 
model B as described in Sect. l4.3l bv sampling the ^-distribution 
around a# and is based on ~1000 simulated SEDs in the vicin- 
ity of this parameter set. We chose f ~0.5 ■ x\ an d deduced 
the asymmetric intervals given in Table [5] The reader is cau- 
tioned to not confound these confidence intervals with a clas- 
sical lcr error interval because the model is highly degenerated 
and consequently we are not dealing with independent parame- 
ters. Additionally, the selection of A;^ f is based on the subjec- 
tive quality of the fit to the observational data as Xsed red * ■ 
Comparison of the confidence intervals with the parameter sub- 
space defined in Table|3]shows that the upper interval of R oa has 
been clipped. As discussed in Sect. 5. 1, we cannot give an upper 
constraint for this parameter with available observations and ac- 
cordingly did not explore the confidence interval beyond 500 AU. 
The parameters r er r and L„ describing the source might also be 
influenced by clipping because their intervals almost cover the 
given subspace. The error bounds reflect the degeneracy in the 
vicinity of our best-fit model B and remind us of the difficulties 
of constraining a model in a high-dimensional parameter space 
using only SED data. 
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Fig. 8. Direct comparison of horizontal cuts through the bright- 
ness map at X = 1.3mm of the deconvolved PdBI observation 
and our best-fit models. 
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Fig. 9. Direct comparison of the synthesized SEDs of our best-fit 
models with observations in the investigated wavelength domain 
from 5/zm to 3mm. 



5.3. Fitting of the SED with a simple flared disk 

The best-fit model C has a slightly steeper and colder radial tem- 
perature distribution than the previous model B with g~0.43 and 
Jioo~12K for distances /?>5AU and therefore lies closer to the 
molecular line observations by P06, see Fig. [7] Although its SED 
appears to be a good match to observation, the 1.3mm bright- 
ness map obtained from model C does not show a central dip, 



see Fig. [8] This finding suggests that it is improbable to simulta- 
neously reproduce the SED and PdBI observations with a simple 
flared disk ansatz. 

5.4. Simultaneous fitting of SED and mm map 

The tradeoff model D is similar to parameter set B, apart from 
its smaller outer radius. The radial temperature distribution in 
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the mid-plane is close to model B with q~0.36 and rioo~19K 
for 7?>5AU. All deduced properties are comparable except for 
its mm brightness map, which features an improved shape and 
peak flux although the contrast of the double peak is still insuf- 
ficient, but the differences are not significant. The optical depth 
of model D at A = 1.3mm is t~1.8 in the mid-plane and t~0.9 
along the line of sight, so the disk is semitransparent in the mm- 
wavelength domain. 

5.5. Discussion 

Our study confirms previous results using an unbiased, extended 
disk model and a completely different minimization method, 
which is less susceptible to secondary minima. Nevertheless, 
none of the Markov chains of run B converged toward rj~l, al- 
though W02 presented a well-matched simple flared disk model. 
The two domains with small and large tj appear to be separated 
by a "ridge" in the x 1 -distribution that our specific implementa- 
tion of the annealing schedule cannot overcome easily within the 
performed number of steps. To check the possibility of a simple 
flared disk solution, we fixed rj- 1 and found model C and similar 
parameter sets with^ ED comparable to model B. Their appear- 
ance at 1.3mm suggests that the deficit of mm emission toward 
the center of HH 30 is unlikely to be explained by an optical 
depth effect in a nearly edge-on simple flared disk, see Fig. [8] 
Therefore a change of mm opacity by a factor of >10 at ~45AU 
between an inner and outer domain is required, as discussed in 
G08. This result can be directly understood from the mm bright- 
ness image. 

The typical dust temperature 7d us t at 50AU in our best-fit 
models is >20K, comparable to temperatures deduced from 
molecular line observations (P06), see Fig. [7] The surface 
brightness 7b of an optically thick emitter filling the whole 
beam is just its mean temperature, therefore we can deduce 
from the actually observed peak brightness of T\, <0.5K that any 
optically thick emitter has to be geometrically diluted by a factor 
of 7dust/?b ^40. The orientation of the beam resolves the disk 
of HH 30 radially and leaves only vertical settling of dust as an 
explanation for an optically thick scenario. The dilution factor is 
accordingly given by the ratio of the projected disk thickness to 
the vertical resolution, 0.59"~83AU at the distance of ~140pc 
for Taurus. Thus, the apparent disk thickness should be ~2AU, 
incompatible with the derived inclination of ~83°, because 
sin7°~0.12 and the diameter of the optically thick parts should 
be in the order of the distance between the peaks, i.e. ~100AU. 
To still be able to explain the deficit of mm emission towards 
HH 30 by self-absorption, two effects have to be combined: 



1 . A very strong settling of larger grains to obtain a scale height 
/i~lAUatfl~50AU. 

2. A substantial disk warp aligning the inner parts to appear 
nearly edge-on (>88°), while the outer disk, dominating the 
scattered light images, remains at an inclination of ~83°. 

While this scenario cannot be excluded, these peculiar fea- 
tures are beyond the scope of models tested so far, because all 
assume cylindrical symmetry. Given that inner holes have now 
been discovered in many disks, e.g. in LkCa 15 (Pietu et al. 
2006), GM Aur (Dutrey et al. 2008; Hughes et al. 2009), and 
CB 26 (Sauter et al. 2009), a low-opacity cavity is a simple 
and viable explanation for the appearance of HH 30 at mm- 
wavelengths. Reanalysis by G08 of the I3 CO data published in 
P06 shows that the lack of emission for >3km/s can be explained 



by removing gas from the inner (33+7)AU. Anglada et al. (2007) 
proposed a wide binary to explain the jet wiggling that would 
clean the inner disk by tidal interaction as exemplified by G08. In 
this study we showed that a good fit of the SED and the double- 
peaked mm map can be achieved simultaneously if a depletion 
zone is included (Model D). Hence, this notion is now supported 
by several independent pieces of evidence. 

Attempting to account for all known broadband measure- 
ments confirms a number of previously noted facts, and helps 
revealing their possible origins. Under our hypothesis of ho- 
mogeneous dust properties, no single model can reproduce the 
ensemble of available data. A major problem is the outer ra- 
dius of the dust disk. Resolved observations in the NIR and 
the SED are consistent with /? out ~200..500AU, while the con- 
tinuum visibilities measured with the PdBI alone clearly indi- 
cate R out = (135+15)AU, a discrepancy that cannot be explained 
with a single dust species and suggests a radial dependency of 
the dust composition. A recently conducted high-resolution sur- 
vey at mm-wave lengths with a sample of 10 TTauri disks found 
the same radial change in opacity (Guilloteau et al. 2011), in- 
dicating a negative gradient of mean grain size d{a) /dR<0 in 
the outer disk with populations of smaller particles present in 
this region. As pointed out by G08, the outer radius of the gas 
disk (~420AU) is consistent with this interpretation. Further op- 
timization should be continued in a follow-up study when addi- 
tional high-resolution observations in the (sub)mm domain are 
available that reveal the spatial dependence of the spectral index 
/3mm- Our results show that the mm observations described in 
G08 can be reproduced by a deficit of emission within ~45AU, 
while the SED indicates in the MIR substantial amounts of warm 
dust in this region. We can explain the lack of flux from the in- 
ner zone at 1.3mm by an overall depletion of dust. However, 
the large grain size required to produce sufficient mm flux is 
inconsistent with chromaticity observed in scattered light im- 
ages. Dust settling may explain this effect, because the scale 
height of large, and therefore mm emitting, grains is only loosely 
constrained. Another equally likely scenario is the migration of 
small dust grains coupled to the gas flow from the outer ring 
through the inner region, e.g. along streamers (Grinin et al. 2004; 
Kley et al. 2008; Tambovtseva et al. 2008), while larger grains 
decouple from this motion and are somehow effectively removed 
in the depletion zone. Either way, the inner region is not devoid 
of gas and dust as the SED in the MIR and the active jet clearly 
demonstrate. Only interferometers with sufficient sensitivity and 
sub-arcsecond resolution like ALMA will be able to significantly 
reduce additional model degeneracies introduced by parameters 
describing an inhomogeneous grain size distribution, and disen- 
tangle whether vertical stratification is required in addition to 
radial segregation. 

6. Summary 

We conducted a parameter study to find an unbiased, self- 
consistent model of HH 30's disk able to predict all available 
continuum data. The investigation started with the c-disk model 
and a standard dust composition as used in previous modeling ef- 
forts. Preliminary results exhibited incongruent fluxes and spec- 
tral indices in the MIR and for (sub)mm-wavelengths, leading to 
a refinement of the model. To compensate for the discrepancies 
between model predictions and observations, we added a vari- 
able grain size distribution and an attenuated region, thereby in- 
troducing larger grains and warm dust near the central source. 
The introduction of a depleted region and variable dust size 
added five parameters to our study and rendered an exhaustive 
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search for the global best fit impractical, which led us to imple- 
ment SA. Modification of the density by an attenuation 77 was 
mainly selected for its simplicity and many configurations com- 
patible with the observed flux are conceivable. Our model con- 
centrates most of the dust near the mid-plane, causing higher op- 
tical depths in this region. Mid-infrared flux is hence mainly pro- 
duced near the inner edge and in warm, upper layers because the 
inner parts toward the mid-plane are optically thick in this wave- 
length domain. Less massive configurations with lower optical 
depth may reproduce the observed flux in the MIR, therefore the 
deduced dust mass of the depleted region should be interpreted 
as an upper limit. 

Earlier work reported extensive ambiguities when fitting ob- 
servations, even to the point that a single best-fit parameter set 
could not be given (W04), and we encountered the same obsta- 
cle, especially when fitting single data sets. 

The parameter set A found during optimization of NIR im- 
ages reproduces the appearance of HH 30 in scattered light, but 
is not able to correctly predict any observation at longer wave- 
lengths. Our inability to reproduce the contrast between sec- 
ondary maximum and the obscuring belt in the NIR images may 
indicate the presence of additional features not present in our 
simple model, such as warps, an upwhirled dust cone, or asym- 
metrical illumination. 

Degeneracies were also present during SED fitting, though 
new data sets have substantially improved the constraints on our 
predictions, most notably the observations by IRS/SST in the 
MIR. The parameter set of model B was found during optimiza- 
tion of the SED and replicates the flux in MIR quite well, but 
does not exhibit the correct appearance in the NIR and its mm 
brightness map is too extended, showing low peak brightness 
and little contrast. The overall agreement of observation and sim- 
ulation at mm- wavelengths is only satisfactory for this model be- 
cause the fluxes at 2.7mm and 3.4mm are underpredicted. We de- 
duced confidence intervals for this model using Markov chains 
that sampled the vicinity of a^, see Table|5] 

To improve the model and reduce degeneracies we com- 
bined the optimization of two different classes of observations, 
photometry from MIR to mm-wavelengths and interferomet- 
ric imaging at 1.3mm, and optimized with our new staggered 
MCMC method that does not weigh individual x\ during opti- 
mization. This approach did remove some degeneracies and the 
resulting best-fit model D has an improved mm brightness map 
without compromising the SED. 

The findings of our parameter study extend the results of pre- 
vious works, see Tables |2] |4] and|5]for comparison. This is not 
self-evident because our modeling effort was not biased by fix- 
ing parameters a priori and the standard ansatz was modified by 
adding an attenuated inner zone, which introduced new degen- 
eracies. 

A greater refinement of the model appears to be necessary 
to satisfactorily explain all observational data at once. The most 
prominent discrepancies between model A and the other best- fit 
models lie in the grain size distribution and hence dust mass and 
chromaticity. This suggests that a disk of larger settled grains 
embedded in an overlaying disk of smaller dust grains may be 
able to reconcile theory and observation. A model accounting 
for spatially inhomogeneous grain size distributions has to in- 
clude at least two different dust species and several new pa- 
rameters, e.g. D'Alessio et al. (2006), Sauter & Wolf (2010). 
Unfortunately, at the moment we lack observational data to dis- 
entangle the additional degeneracy, which significantly reduces 
the plausibility of any such model. Nevertheless, our results are 
consistent with an inner depletion zone of (45+5)AU radius, 



slightly larger than the hole proposed in G08, and suggest dust 
growth, vertical settling, and radial segregation as mm opacity 
changes significantly at ~140AU. 

High-resolution imaging of HH 30 in different wavelength 
domains with next-generation observatories, such as the interfer- 
ometer ALMA and the space-born JWST, will break the inherent 
model degeneracy and help to uncover the structure of the inner 
region, map the spatial dependency of dust properties, and refine 
our understanding of this fascinating object and its evolutionary 
stage in the coming years. 
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